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Abstract 

We study minimal time strategies for the treatment of pollution of large volumes, such as lakes 
or natural reservoirs, with the help of an autonomous bioreactor. The control consists in feeding the 
bioreactor from the resource, the clean output returning to the resource with the same flow rate. We 
first characterize the optimal policies among constant and feedback controls, under the assumption of a 
uniform concentration in the resource. In a second part, we study the influence of an inhomogeneity in 
the resource, considering two measurements points. With the help of the Maximum Principle, we show 
that the optimal control law is non-monotonic and terminates with a constant phase, contrary to the 
homogeneous case for which the optimal flow rate is decreasing with time. This study allows the decision 
makers to identify situations for which the benefit of using non-constant flow rates is significant. 
Keywords. Environmental engineering, biotechnology, waste treatment, continuous systems, minimum- 
time control. 

1 Introduction 

The fight against eutrophication of lakes and natural reservoirs (excessive development of phytoplankton 
associated with an excess of nutrients) constitutes a major challenge. Such an ecological question has given 
rise to many studies over the last 30 years (see, for instance, the surveys [5] or [25] and references herein). 
To remediate to eutrophication, many techniques such as bio-manipulation or ecological control have been 
proposed with mitigated results. A common point of the proposed remediation approaches is that they are 
usually based on " biotic" actions on the lake trophic chain dedicated to the restoration of the equilibrium 
of the local ecosystems. To do so, most studies are based on empirical knowledge. However, since the 
seventies, the use of eutrophication models (from heuristic data-based models at steady state to more recent 
dynamical mass-balance based models) together with optimal control techniques have been proposed (cf. [B] 
and references herein). 

In the present paper, an alternative to these techniques is studied using a very simple model of the 
lake. It is assumed that a small bioreactor is available to treat the polluted water in removing a substrate 
considered as being in excess in the lake water. More particularly, we consider a natural water resource 
of volume V polluted with a substrate of concentration Si. As underlined above, typical examples of such 
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natural water resources to be treated are lakes or water tables that have been contaminated with diffused 
pollutant as organic matter or nutrients. The objective of the treatment is to make the concentration of 
such pollutant/contaminant Si decreasing down, as fast as possible, to a prescribed value Sj, with the help 
of a continuous stirred bioreactor of volume V r . The reactor is fed from the resource with a flow rate Q, and 
its output returns to the resource with the same flow rate Q, after separation of biomass and substrate in 
a settler (see Figure [I| . The settler avoids the presence of excessive biomass used for the treatment in the 
natural resource, that could bring undesirable sludge and possibly lead to an increase of the eutrophication. 
We assume that during the whole treatment, the volume V of the resource does not change. 

Q 

I ~~1 

bioreactor 
- V r 
Q 

I settler 
biomass 

Figure 1: Interconnection of the bioreactor with the resource. 

Since the pioneer work by 5], the optimization of bioreactors operation has received a great attention in 
the literature; see [19j [3j [2] for reviews of the different optimization techniques that have been used in 
bioprocesses. Among them, the theory of optimal control has proved to be a generic tool for deriving 
practical optimal rules HS1 HI] • Clearly, one can distinguish two different kinds of problems depending 
on the continuous or discontinuous operation mode of the process. On one hand, if the process is operated in 
fed-batch, the control objective is usually to optimize trajectories for attaining a prescribed target in finite 
time or maximizing the production at a given time 0[Tg[T2[TTJ[T![]ni[2S[T^ On the other 

hand, the optimal control of continuous processes usually involves a two steps procedure. First, the optimal 
steady state is determined as a nominal set point, maximizing a criterion [271 126) . The benefit of operating 
a periodic control about the nominal point can be analyzed [TJ HD] . Then, a control strategy that drives the 
state about the nominal set point from any initial condition is searched for [13., possibly in the presence of 
uncertainty on the model using extremum seeking techniques [301 EH EE |4] 

Concerning these strategies, the problem studied in the present paper exhibits several original points 
with respect to the contributions available in the literature. Indeed: 

- The actual control problem is dedicated to the optimization of transient trajectories - as in the case of 
fed-batch processes - while it is actually a continuous process. It is due to the fact that in a standard 
optimal minimal-time problem of a bioprocess, the volume of water to be processed is completely 
decoupled from the bioreactor. In other terms, the problem is to process, using a biological reactor, a 
given volume of " substrate" which is finally released in the environment after processing (whatever it 
is operated continuously or discontinuously) . In the present problem, the treated water is immediately 
recycled into the lake. From the modelling point of view, this introduces an original coupling via the 
dilution of the treated water with the polluted one. 




2 



- The lake and the reactor are isolated in the sense no biomass is supposed to be present in the water 
resource. The biomass used as a catalyst in the bioreactor is separated from the treated water and 
withdrawn from the overall process. Thus, in particular, the quantity of available biomass is not a 
limiting parameter. 

We consider the usual chemostat model for describing the dynamics of the bioreactor: 

Sr = —/J,(Sr)X r + —(Sl — Sr) 

• Q W 

X r — /i(S r )X r — jrX r 

where S r and X r stand for the concentrations of substrate and biomass, respectively. For sake of simplicity, 
we assume that the yield coefficient of this reaction is equal to one (at the price of changing the unitary value 
of the biomass concentration, that is always possible). The growth rate function /i(-) fulfills the properties 

Assumption Al. 

a. Function /x(-) is increasing and such that fi(0) = 0. 

b. Function /i(-) is concave. 

A reasonable hypothesis is to assume that the volume of the resource is much larger than the bioreactor 
one: V >> V r , and that the possible variations of the manipulated variable Q are slow compared to the 
time scale of the bioreactor dynamics. Consequently, one can consider that dynamics ((T|) is fast and its 
trajectories at the quasi-steady state (S*,X*) — (S r (Q), Si — S r (Q)), where S r (Q) fulfills fi(S r (Q)) — Q/V r 
(see the usual equilibria analysis of the chemostat |24| ) . 

Problem: The optimization problem consists in driving in minimal time the concentration of the resource 
down to a prescribed value S_ t > 0, playing with the control variable Q > 0. In Section 2, we assume that this 
concentration is uniform in the resource, while in Section 3 we study the effect of a spatial inhomogeneity. For 
each case, we characterize the optimal policy Q* (resp. Q opt {-)) among constant (resp. feedback controls). 
Section 4 is devoted to numerical simulations and discussions. 

2 The homogeneous case 

The dynamics of the resource concentration is simply 

St = ~(S r (Q)-St). (2) 
Notice that under Assumption Al.a, choosing Q is equivalent to choosing S r as a control variable: 

Si = an(S r )(S r - Si), S r e (0, Si) (3) 

where we denote a — V r /V. 

Proposition 1 Under Assumption Al, the best constant control Q* is defined as Q* — V r fi(S*), where S* 
is the unique minimum of the function 

on the interval (0,S_i). 
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Proof. For a constant control, solutions of can be made explicit: 



Si(t) = S r (Q) + (Si(0) - S r {Q))e-^ 



(5) 



as well as the time Tf(S r ), given in (UJ), for reaching the target with Q = V r fi(S r ). The function !/(•) tends 
toward +00 when S r tends toward or S_ t . Consequently, its infimum is reached on the interval (0,5;). 
Denote by Tj its minimum, that we fix in the following. Then, for each constant control S ri one has 



dSi(Tf) 



dS r 



i + a/i'(^)r;(Sz(o)-s r ) 



-oi*(s r )r; 



d?Si(T*) 
dS? 



2^'(S r ) + (afi'(S r ) 2 T} - M "(5 r )(5,(0) - S* r ) 



aT * e - a »(s r )T} 



and one deduces with Assumption Al that the map S r n> Si(Tj) is strictly convex. Notice that one 
has necessarily Si(Tf) > S_i and Si(Tj) = S_ t when S r = S* realizes the minimum of the function Tf(-) 
Consequently, the optimal control S* is unique. 



□ 



Proposition 2 Under Assumption Al, the optimal feedback fulfills Q opt (Si) — V r [i(S° pt (Si)) with 

S^iSt) G argmax/i(S,.)(S/ - S r ) . 

s r e(o,Si) 



(6) 



Moreover, t 1— > Q opt (t) is decreasing along any optimal trajectory. 

Proof. It is straightforward to check that the optimal feedback S° pt is the one that makes the time deriva- 
tive of Si, given by ((3]), the most negative at any time. A necessary condition is to have have fi'(S° pt )(Si — 
S° pt ) = ii{S° pt ). Deriving this last expression w.r.t. time, one has S° pt (2fi' (S° pt ) + fi" (S° pt )(S° pt - S t )) = 
n'(S? pt )Si and from Assumption Al, on obtains S° pt < 0. □ 

For usual growth functions, one obtains the expressions: 



linear: /i(s) = /is, 


Monod: 


M ( S ) = 


K + s' 


S° pt {Si) = Si/2 


S° pt (Si) = 


Vk 2 - 


YKSi-K 



3 Consideration of a spatial inhomogeneity 

The simplest way to introduce a gradient of concentration in the model of the resource is to consider two 
compartments of volumes V\ , V2 such that V = V\ + V2 (see Figure [5]) , that we assume to be large with 
respect to V r . Water is pumped from the first one while the clean one is rejected in the second one. 



Denoting cti — V r /Vi (i — 1, 2), one obtains the dynamics 

51 = - Si) = am(S r )(S2 - Si) 

1 (7) 

Q 

5 2 = — (S r (Q) - S 2 ) = a 2 /J,(S r )(S r - S 2 ) 

v 2 

and can easily check that the domain V = {(Si,S 2 ) £ \Si > S 2 } is invariant for any control S r (-) 
such that S r (t) S (0,S 2 (t)} for any t > 0. We shall consider initial conditions in T> and define the target 
T={(Si,S 2 )eV\Si<S l }. 
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bioreactor 



settler 



biomass 



Figure 2: Consideration of non- homogeneity in the resource. 



For p G [0, 1] and r > 0, we define the function 



A(p, r) = 



(1 — p)e 1 -p 



pe 



1 - 2p 



and for Sq > Si the function 



B(S r ) = 



1+ 2 
A(l-p,r) 



if p = 

ifpe (0,|) 

ifp = i 
ifpe (ii] 



S r g (0,5,) 



5o — 5 r 

Proposition 3 Let p = Vi/V. For initial conditions such that 5i(0) 



52(0) = 5q > S_i, the best constant 



control Q* and time to reach the target T are defined by Q* — V r fi(S*), where S* is such that the graph 
of B(-) touches tangentially the graph of S r 1— > A(p, ix{S r )TJ) at S r = S*. 

Proof. The solution of (0 with constant control S r £ (0, 5;) can be made explicit: 

S 1 (t) = S r + (S -S r )A(p,fi(S r )t) , 

and the time X/ to reach the target fulfills A(p, fi(S r )Tf) = B(S r ). Notice that one has the property 

Sl(t) >S t ^ A(p, fl(Sr)t) > B(S r ) , 

from which the statement of the proposition follows. □ 

If one consider the family of functions AT(S r ) = A(p, (i(S r )T), parametrized by T > 0, one has a graphical 
interpretation of the optimum depicted in Figure [3] 

£?(•) is a concave function and one can check that At(-) are strictly convex for values of S r large enough 
(under Assumption Al). Consequently, there cannot exist more than one best constant control S* in the 
domain where At* is convex. 

Remark 1 Notice that the case of an homogeneous resource can be seen formally as the limiting case p = 
( although cases when p is close to are not compatible with assuming that V\ , V2 are large with respect to 
V r ). One can easily check that A(p, r) < A(0, r) for sufficiently large values of t. Furthermore, times 
are increasing with respect to Sq — 5,. Consequently, for initial conditions such that Sq is far from 5,, the 
time Tj is larger for an homogeneous resource than for non-homogeneous one, even when the parameter p 
is unknown and the control S* is determined for the homogeneous case. 
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Figure 3: Graphical determination of S* and TJ. 
For 5i > 52 > and 7 > 0, we define 



(f>(Si,S 2 ,^,S r ) = fi(S r ) 



1 + l (s 2 -s r ) 



{Si - s 2 ) 



(8) 



V>(Si,S 2 , 7 ) = ^)"7^. (9) 
The proof of the following lemma is left to the reader. 

Lemma 1 Under Assumption Al, for S\ > S 2 > and 7 > 0, the function (t>(Si, S 2l 7, ■) is strictly concave 
on [0, S2] and the property 

max (f>(Si,S 2 ,7,S r ) = </>(5i,5 2 ,7,5 2 ) 
s r e(o,s 2 ] 

is fulfilled exactly when tp(S±, S 2 , 7) > 0. 

Proposition 4 Under Assumption Al, from any initial condition in T> \ T, the optimal control Q opt (-) 
consists in reaching a subset X C T> \ T from which the constant control Q opt — V r fJ,(S_ 2 ) is optimal until 
Si(-) reaches Sj, where S_ 2 is the value of S 2 when X is reached. Moreover, t i-> Q opt (t) is increasing when 
approaching the set X. 

Proof. Recall first that T> is invariant. If S\ = S 2 > S_ii the feedback S r — S 2 cannot be optimal (this 
would imply Si = S 2 = at any time). So, any optimal trajectory is such that Si(t) > S*2(t) for any t > 0. 
Let us write the Hamiltonian, along with the adjoint equations: 

H = — 1 + max n(S r ) [ai\i(S 2 - Si) + a 2 \ 2 (S r - S 2 )} 
s r e[o : s 2 ] 

A a = a lf i(S°P t )X 1 , \i(T opt ) < 

A 2 = ^(S? pt ){a 2 X2 - aiAi) , X 2 {T opt ) = 

One deduce immediately that Ai (t) < for any t > and can consider the function 

lit) = (10) 
aiXi{t) 

that fulfills 7 = ii(S° pt ) [(a 2 - «i) 7 - a 2 ] , j(T opt ) = 0. 

Notice that 7 = implies 7 < and then one obtains 7(t) > for any t £ [0, T opt ). 



(i 



When Si > S2, optimizing the Hamiltonian is equivalent to maximizing <p(S±, S2, 7, •) (defined in ©), 
and then Lemma Q] provides the uniqueness of S° pt . A straightforward calculus gives 
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dt \Si — S2 



Si — S2 



1 S2-S?* \ 1 
Si — S2 ) 



From j(T opt ) = we deduce the existence of i 6 [0, T opt ) such that ^ \ ~ST 
for ip given by ([9]), one has 



ip = S 2 [ f/'(S 2 ) - 



Si 



7 



y Si — S2 



< for all t S [i, T„pt]. Then, 



MOS2), 



which is positive for t G [t, T op t]. Since ip > at T op t, there exists i s G [t,T op t) such that ^ > for t > t s . 
Defining t = inf{i s G [i, T op {) : ip > for i G [i s ,T pt]}, with Lemma [1] one concludes that the optimal S° pt 
is constant equal to S_ 2 = Saft) at any time t > i. 

When f > 0, let us write S° pt = US2 with u G [0,1]. The left derivative u{t~) has to be positive and 
S 2 (t) = 0. This implies to have S° pt {t-) > 0. □ 

Proposition 5 Under Assumption Al, for any initial condition inT>\T , the optimal trajectory is unique. 

Proof. We recall, from the proof of Proposition!!! that along any optimal trajectory, one has Si(t) > S2 (t) 
and 7(f) > for any t G (0, T/). Then, one has Si < and can re-parametrize the dynamics of variables S2 
and 7, defined in (|10[) , in terms of Si instead of time t, and write the non-autonomous dynamics for optimal 
trajectories: 

dS 2 _ a 2 {S° p - S 2 ) . _ . 

-77T — TTf FTT - ' ^2(04) — 02{J- opt) 

dSi ai {S2-S1) 

dj (Q!2 - Q!l)7 - «2 , c n n 

-77T = TTf o~\ 1 7(j2i) = 

where S 1 ^* is the unique maximum of <p(Si, S2, 7, •) on [0, S2]. When S*° pt < S2 and 7 > 0, one has 

f i'(S° pt ) = (12) 



dS r 



S 2 - gg \ 

S*l — S2 Si — S2J 



and hence 



c9 2 



m(S7 4 ) 



< 0. Fix Si, S2 and consider S° pt as a function of 7. From (fT2")l . one has 



d 2 



d 2 6 dS opt 



dS r dj dS 2 r dj 



i)S" 



and from the strict concavity of <fi(Si, S2, 7, •)) given by Lemma[TJ one deduces g ^ 
The Jacobian matrix of system (lll[) is of the form 



< 0. 



Ct2 



ggapt 



((ct2-Qi)7-"2) 
ai(S 2 -Si) 2 



ai(Sa-Si) 97 



from which one observes the non-negativity of off-diagonal terms, because («2 — c*i)7 — &2 — 7 / ' l i {S° pt ) < 0. 
So, the dynamics (fTTj) is cooperative (in time Si), and since 7(5;) = 0, one deduces that two solutions of 
(fTTjl cannot cross in the (Si, 62) plane. Finally, one obtains the uniqueness of the optimal trajectory for a 
given initial condition in T> \ T. □ 
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Remark 2 When Si(0) = S 2 (0) = Si, from the expression of the Hamiltonian we obtain that the optimal 
control S° pt is such that at the beginning it maximizes fi(S r )(Si — S r ). Therefore, it is exactly the same as 
in the homogeneous case of Proposition [H Measuring the initial rate of variation of S 2 gives an estimation 
of the parameter a 2 to fit the model, as one has 

S 2 (0) = a 2 MS r opt )(S 2 (0) - S° r pt ) . 



If it is close to a, then the model with one compartment should suit. 
For Si > Sj > S 2 , we define when ai 7^ a 2 : 



f (Q° Q°\ — V Vgj " " ' " ' R — I 1 

[a 2 - oli)\Pi - S^) \a 2 



and when a, = a 2 : / (S?, S 2 °) = -^fE^l, = e . 

Proposition 6 The set I, where a constant control is optimal, is given by 

1 = {(S°i, S° 2 ) e (g„+cx>) x (0,SJ s.t S 2 ° < S 2 or fi(S°)fo(S°i, S° 2 ) < ,/(S 2 °) } 
where S2 is the unique solution in (0,Sj) of 

»(S 2 )=pLi , (S 2 )(St-S 2 ) . (13) 



Proof. With control Q — y r /i(S 2 ), S 2 (-) is equal to S 2 and solution Si(-) can be made explicit. Then, 
time Tf such that Si(Tj) — S z , and solution 7(-) such that j(Tf) = can be also made explicit. Let 
f(t) = j(t)/(Sx(t) — 5*2). According to Lemma[TJ this constant strategy is optimal exactly when 

//(S°) > (i(S$)f(t) , te[Q,T f ]. (14) 

One can easily check that / = a 2 /x(S 2 )(7 — l)/(Si — S 2 ), and consequently, / cannot be null more than one 
time (recall from the proof of Proposition [4] that 7 is non- increasing). One has also /(0) > 0, f(Tf) = 0, 
and f'(Tf) < 0. 

If /'(0) < 0, then condition (|14[) is equivalent to have (5°, S 2 ) below the graph of the curve C defined by 
n{S%)fo{s%,S%) = (J.'(S$). A straightforward but lengthy computation gives /(0) = /o(S?, S%)). One can 
also check that /'(0) < is equivalent to have (S^S^below the line C defined by S? - S% = /3(S, - S$). 
The intersection point (Si, S 2 ) of C and £ is given by S 2 solution of (fl"3]l . its uniqueness being guaranteed 
by the concavity of fj,. One can easily check that C is below C for any S 2 6 [S 2 , S z ]. 

If /'(0) > 0, on can check that maxtf(t) = 1//3(S/ — S 2 ) and then condition (fT4| is equivalent to have 
S 2 < S 2 . Moreover, the straight line S 2 = S 2 is below the graph of the curve C in the interval [S ; , Si] (see 
Figure gj. □ 

For the Monod law, one can find 

S 2 = \ (-K(l + /3) + VK 2 (1 + P) 2 + 4K/3S ; ) . 



4 Discussion and numerical simulations 

The benefits of our theoretical analysis is to identify efficient pumping strategies and some of their robustness 
properties. We summarize those contributions in terms of the following rules for the decision makers: 
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Sf Si S i 

Figure 4: Backward integration of the extremals. 

1. The profit of using the optimal feedback strategy compared to the best constant one, can be easily de- 
termined numerically (see the simulations below). As expected, the more the resource is initially polluted, 
the more the improvement of the feedback policy is significant. Depending on the ratio "initial pollution 
over desired maximal pollution level" , the decision maker can then decide whether it worth adopting a time- 
varying strategy. 

2. A spatial inhomogencity of the pollution concentration improves the treatment time on the condition 
that the resource is enough polluted. Moreover, applying the best constant strategy as if the resource was 
perfectly homogeneous is robust with respect to uncertainty on the inhomogeneity parameter in the sense 
that it provides a guaranteed time (see Remark [1}. 

3. Measuring the initial speed of variation of concentration at two remote locations in the resource allows 
to identify the inhomogeneity parameter of the model (see Remark [2]). Then, the decision maker can decide 
if it worth considering a feedback strategy with two measurement points instead of one. 

4. The optimal feedback strategy for the inhomogeneous case consists in applying a constant flow rate when 
the concentration S2 reaches a prescribed value given by Proposition [SJ The concentration S% is then main- 
tained constant, without having to measure the concentration S± (see Figure 

Simulations have been conducted for the Monod law with /i max = Is , K = lmol.m~ 3 , and volumes 
V = 1000 m 3 , V r — lm 3 . The initial concentration of pollutant has been chosen equal to lmol.m~ 3 
uniformly. Figure [5] shows the comparison of minimal times for different values of S_ t (the curves corresponds 
to different values of the parameter p) . 

On this example, one can see that for S_[ = 0.01, the minimal time among constant controls is about twice 
larger than among feedbacks. The influence of inhomogeneity is also quite significant. 

The optimal feedback ([6]) for the homogeneous case is a simple law that provides a decreasing flow rate 
Q w.r.t. time (cf. Proposition [2]) , contrary to the inhomogeneous case for which it is non-monotonic (cf. 
Proposition|4|). In Figure [6l we compare the optimal policy Q opt (-) for p — 0.4 and S_i — 0.1 with Qi(-), 
resp. Qi{-) applying the formula ([6} on measurement Si, resp. S%. 

The true optimal feedback control, in the model that consider both measurements, is more sophisticated 
in the sense that it anticipates the approach to the target, increasing the flow rate and freezing it. 
The study has been made assuming that the steady-state characteristics Q i-> S r (Q) of the bioreactor is 
perfectly known. Uncertainty on this map as well as on measurements will be the matter of a forthcoming 
work. 
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Figure 5: Comparison of optimal times. 
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Figure 6: Optimal and sub-optimal policies. 
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